Here I study how biomass density changes across treatments in the PatchSizePilot. In particular, I’m studying how biomass density changes across:
culture_info = read.csv(here("data", "PatchSizePilot_culture_info.csv"), header = TRUE)
load(here("data", "population", "t0.RData")); t0 = pop_output
load(here("data", "population", "t1.RData")); t1 = pop_output
load(here("data", "population", "t2.RData")); t2 = pop_output
load(here("data", "population", "t3.RData")); t3 = pop_output
load(here("data", "population", "t4.RData")); t4 = pop_output
load(here("data", "population", "t5.RData")); t5 = pop_output
load(here("data", "population", "t6.RData")); t6 = pop_output
load(here("data", "population", "t7.RData")); t7 = pop_output
rm(pop_output)
#Column: time
t0$time = NA
t1$time = NA
#Column: replicate_video
t0$replicate_video = 1:12 #In t1 I took 12 videos of a single
t1$replicate_video = 1 #In t1 I took only 1 video/culture
t2$replicate_video = 1 #In t2 I took only 1 video/culture
t3$replicate_video = 1 #In t3 I took only 1 video/culture
t4$replicate_video = 1 #In t4 I took only 1 video/culture
t5$replicate_video = 1 #In t5 I took only 1 video/culture
t6 = t6 %>%
rename(replicate_video = replicate)
t7 = t7 %>%
rename(replicate_video = replicate)
#Create an elongated version of t0 so that each of the 110 cultures can have 12 video replicates at t0.
elongating_t0 = NULL
for (video in 1:nrow(t0)){
for (ID in 1:nrow(culture_info)) {
elongating_t0 = rbind(elongating_t0, t0[video,])
}
}
ID_vector = rep(1:nrow(culture_info),
times = nrow(t0))
elongating_t0$culture_ID = ID_vector
#Merge previous data-sets
t0 = merge(culture_info,elongating_t0, by="culture_ID")
t1 = merge(culture_info,t1, by = "culture_ID")
t2 = merge(culture_info,t2, by = "culture_ID")
t3 = merge(culture_info,t3, by = "culture_ID")
t4 = merge(culture_info,t4, by = "culture_ID")
t5 = merge(culture_info,t5, by = "culture_ID")
t6 = merge(culture_info,t6, by = "culture_ID")
t7 = merge(culture_info,t7, by = "culture_ID")
ds = rbind(t0, t1, t2, t3, t4, t5, t6, t7)
rm(elongating_t0, t0, t1, t2, t3, t4, t5, t6, t7)
ds$time_point[ds$time_point=="t0"] = 0
ds$time_point[ds$time_point=="t1"] = 1
ds$time_point[ds$time_point=="t2"] = 2
ds$time_point[ds$time_point=="t3"] = 3
ds$time_point[ds$time_point=="t4"] = 4
ds$time_point[ds$time_point=="t5"] = 5
ds$time_point[ds$time_point=="t6"] = 6
ds$time_point[ds$time_point=="t7"] = 7
ds$time_point = as.character(ds$time_point)
ds$day = NA
ds$day[ds$time_point== 0] = 0
ds$day[ds$time_point== 1] = 4
ds$day[ds$time_point== 2] = 8
ds$day[ds$time_point== 3] = 12
ds$day[ds$time_point== 4] = 16
ds$day[ds$time_point== 5] = 20
ds$day[ds$time_point== 6] = 24
ds$day[ds$time_point== 7] = 28
#Column: eco_metaeco_type
ds$eco_metaeco_type = factor(ds$eco_metaeco_type,
levels = c('S',
'S (S_S)',
'S (S_L)',
'M',
'M (M_M)',
'L',
'L (L_L)',
'L (S_L)'))
ecosystems_to_take_off = 60 #Culture number 60 because it was spilled
ds = ds %>%
filter(! culture_ID %in% ecosystems_to_take_off)
ds_for_evaporation = ds
ds = ds %>%
select(culture_ID,
patch_size,
disturbance,
metaecosystem_type,
bioarea_per_volume,
replicate_video,
time_point,
day,
metaecosystem,
system_nr,
eco_metaeco_type)
ds = ds[, c("culture_ID",
"system_nr",
"disturbance",
"time_point",
"day",
"patch_size",
"metaecosystem",
"metaecosystem_type",
"eco_metaeco_type",
"replicate_video",
"bioarea_per_volume")]
ds_regional = ds %>%
filter(metaecosystem == "yes") %>%
group_by(culture_ID,
system_nr,
disturbance,
day,
time_point,
patch_size,
metaecosystem_type) %>%
summarise(patch_mean_bioarea_across_videos = mean(bioarea_per_volume)) %>%
group_by(system_nr, disturbance, day, time_point, metaecosystem_type) %>%
summarise(regional_mean_bioarea = mean(patch_mean_bioarea_across_videos))
metaecosystems_to_take_off = 40 #System 40 was the system of culture 60 that I spilled
ds_regional = ds_regional %>%
filter(! system_nr %in% metaecosystems_to_take_off)
We want to understand if the regional biomass produced by an ecosystem with a small and a large patch (metaecosystem_type = S_L) is lower than the regional biomass produced by an ecosystem with two medium patches (metaecosystem_type = M_M).
Let’s start by plotting the single ecosystems to see that everything is fine. To make the patterns clear let’s plot the low disturbance and high disturbance in two different figures.
Let’s then show these data as boxplots.
We can see that the regional biomass was higher for meta-ecosystems with the same size, regardless of disturbance. As both disturbance levels showed the same pattern, I would just keep the one with low disturbance in the publication, as it shows a stronger pattern.
First of all, let’s modify the data set including the regional biomass of our meta-ecosystems. In this data set, we want to have the regional biomass of the meta-ecosystems (averaged first across videos and then across patches) in which we:
Include only the meta-ecosystems in which patches had both medium size (metaecosystem_type = M_M) and meta-ecosystems in which patches had one a small size and the other large size (metaecosystem_type = S_L).
Take off the first two point (day 0 and day = 4). This is because the first perturbation happened only at day 5.
ds_regional_MM_SL_t2t7 = ds_regional %>%
filter (metaecosystem_type == "M_M" | metaecosystem_type == "S_L",
time_point >= 2)
Let’s start from the largest mixed effect model makes sense to construct. This will include:
We are actually going to construct two full models. One with the correlated and one with the uncorrelated random slopes and intercept, the other without
(1) Should we keep the correlation between intercept and slopes (|)?
full_model_correlated = lmer(regional_mean_bioarea ~
metaecosystem_type +
disturbance +
metaecosystem_type * disturbance +
(metaecosystem_type | day) +
(disturbance | day) +
(metaecosystem_type*disturbance | day) +
(1 | system_nr) ,
data = ds_regional_MM_SL_t2t7,
REML = FALSE)
full_model_uncorrelated = lmer(regional_mean_bioarea ~
metaecosystem_type +
disturbance +
metaecosystem_type * disturbance +
(metaecosystem_type || day) +
(disturbance || day) +
(metaecosystem_type*disturbance || day) +
(1 | system_nr) ,
data = ds_regional_MM_SL_t2t7,
REML = FALSE)
anova(full_model_uncorrelated, full_model_correlated)
## Data: ds_regional_MM_SL_t2t7
## Models:
## full_model_correlated: regional_mean_bioarea ~ metaecosystem_type + disturbance + metaecosystem_type * disturbance + (metaecosystem_type | day) + (disturbance | day) + (metaecosystem_type * disturbance | day) + (1 | system_nr)
## full_model_uncorrelated: regional_mean_bioarea ~ metaecosystem_type + disturbance + metaecosystem_type * disturbance + ((1 | day) + (0 + metaecosystem_type | day)) + ((1 | day) + (0 + disturbance | day)) + ((1 | day) + (0 + metaecosystem_type | day) + (0 + disturbance | day) + (0 + metaecosystem_type:disturbance | day)) + (1 | system_nr)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## full_model_correlated 22 1881.2 1942.5 -918.58 1837.2
## full_model_uncorrelated 31 1899.3 1985.7 -918.63 1837.3 0 9 1
Yes (according to BIC and AIC). Interestingly enough, the p-value is 1 (how should I interpret this?).
(2) Should we keep M * D ?
no_interaction_model = lmer(regional_mean_bioarea ~
metaecosystem_type +
disturbance +
(metaecosystem_type | day) +
(disturbance | day) +
(1 | system_nr),
data = ds_regional_MM_SL_t2t7,
REML = FALSE)
anova(full_model_correlated, no_interaction_model)
## Data: ds_regional_MM_SL_t2t7
## Models:
## no_interaction_model: regional_mean_bioarea ~ metaecosystem_type + disturbance + (metaecosystem_type | day) + (disturbance | day) + (1 | system_nr)
## full_model_correlated: regional_mean_bioarea ~ metaecosystem_type + disturbance + metaecosystem_type * disturbance + (metaecosystem_type | day) + (disturbance | day) + (metaecosystem_type * disturbance | day) + (1 | system_nr)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_interaction_model 11 1861.7 1892.4 -919.86 1839.7
## full_model_correlated 22 1881.2 1942.5 -918.58 1837.2 2.5592 11 0.9954
No (according to AIC and BIC). The p-value is also here almost 1 (how should I interpret this?).
(3) Should we keep (M | day)?
no_metaeco_slopes_model = lmer(regional_mean_bioarea ~
metaecosystem_type +
disturbance +
(disturbance | day) +
(1 | system_nr),
data = ds_regional_MM_SL_t2t7,
REML = FALSE)
anova(no_interaction_model, no_metaeco_slopes_model)
## Data: ds_regional_MM_SL_t2t7
## Models:
## no_metaeco_slopes_model: regional_mean_bioarea ~ metaecosystem_type + disturbance + (disturbance | day) + (1 | system_nr)
## no_interaction_model: regional_mean_bioarea ~ metaecosystem_type + disturbance + (metaecosystem_type | day) + (disturbance | day) + (1 | system_nr)
## npar AIC BIC logLik deviance Chisq Df
## no_metaeco_slopes_model 8 1871.0 1893.2 -927.48 1855.0
## no_interaction_model 11 1861.7 1892.4 -919.86 1839.7 15.239 3
## Pr(>Chisq)
## no_metaeco_slopes_model
## no_interaction_model 0.001623 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Yes (according to AIC and BIC - however, BIC is only slightly lower). The p-value is really low (0.002).
(4) Should we keep (D | day)?
no_disturbance_slopes_model = lmer(regional_mean_bioarea ~
metaecosystem_type +
disturbance +
(metaecosystem_type | day) +
(1 | system_nr),
data = ds_regional_MM_SL_t2t7,
REML = FALSE)
anova(no_interaction_model, no_disturbance_slopes_model)
## Data: ds_regional_MM_SL_t2t7
## Models:
## no_disturbance_slopes_model: regional_mean_bioarea ~ metaecosystem_type + disturbance + (metaecosystem_type | day) + (1 | system_nr)
## no_interaction_model: regional_mean_bioarea ~ metaecosystem_type + disturbance + (metaecosystem_type | day) + (disturbance | day) + (1 | system_nr)
## npar AIC BIC logLik deviance Chisq Df
## no_disturbance_slopes_model 8 1860.6 1882.9 -922.31 1844.6
## no_interaction_model 11 1861.7 1892.4 -919.86 1839.7 4.9107 3
## Pr(>Chisq)
## no_disturbance_slopes_model
## no_interaction_model 0.1785
No (according to AIC and BIC - however, AIC is only slightly lower). The p-value is low (0.1785).
(5) Should we keep (1 | system_nr)?
no_system_nr = lmer(regional_mean_bioarea ~
metaecosystem_type +
disturbance +
(metaecosystem_type | day),
data = ds_regional_MM_SL_t2t7,
REML = FALSE)
anova(no_disturbance_slopes_model, no_system_nr)
## Data: ds_regional_MM_SL_t2t7
## Models:
## no_system_nr: regional_mean_bioarea ~ metaecosystem_type + disturbance + (metaecosystem_type | day)
## no_disturbance_slopes_model: regional_mean_bioarea ~ metaecosystem_type + disturbance + (metaecosystem_type | day) + (1 | system_nr)
## npar AIC BIC logLik deviance Chisq Df
## no_system_nr 7 1861.5 1881.0 -923.75 1847.5
## no_disturbance_slopes_model 8 1860.6 1882.9 -922.31 1844.6 2.881 1
## Pr(>Chisq)
## no_system_nr
## no_disturbance_slopes_model 0.08963 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Yes according to AIC (but no according to BIC). The p-value is low (0.09).
Should we keep these parameters?
| Parameter | AIC | BIC |
|---|---|---|
| | (correlation) | Yes | Yes |
| M * D | No | No |
| (M | t) | Yes | Yes (but slightly lower) |
| (D | t) | No (but slightly lower) | No |
| (1 | ID) | Yes | No |
According to AIC, the best model should be M + D + (1 | system_nr). BIC would take out (1 | system_nr). But as we are going for AIC here, our best model should include it.
To make sure everything is okay, let’s take a look at all our models.
anova(full_model_correlated, full_model_uncorrelated, no_interaction_model, no_metaeco_slopes_model, no_disturbance_slopes_model, no_system_nr)
## Data: ds_regional_MM_SL_t2t7
## Models:
## no_system_nr: regional_mean_bioarea ~ metaecosystem_type + disturbance + (metaecosystem_type | day)
## no_metaeco_slopes_model: regional_mean_bioarea ~ metaecosystem_type + disturbance + (disturbance | day) + (1 | system_nr)
## no_disturbance_slopes_model: regional_mean_bioarea ~ metaecosystem_type + disturbance + (metaecosystem_type | day) + (1 | system_nr)
## no_interaction_model: regional_mean_bioarea ~ metaecosystem_type + disturbance + (metaecosystem_type | day) + (disturbance | day) + (1 | system_nr)
## full_model_correlated: regional_mean_bioarea ~ metaecosystem_type + disturbance + metaecosystem_type * disturbance + (metaecosystem_type | day) + (disturbance | day) + (metaecosystem_type * disturbance | day) + (1 | system_nr)
## full_model_uncorrelated: regional_mean_bioarea ~ metaecosystem_type + disturbance + metaecosystem_type * disturbance + ((1 | day) + (0 + metaecosystem_type | day)) + ((1 | day) + (0 + disturbance | day)) + ((1 | day) + (0 + metaecosystem_type | day) + (0 + disturbance | day) + (0 + metaecosystem_type:disturbance | day)) + (1 | system_nr)
## npar AIC BIC logLik deviance Chisq Df
## no_system_nr 7 1861.5 1881.0 -923.75 1847.5
## no_metaeco_slopes_model 8 1871.0 1893.2 -927.48 1855.0 0.0000 1
## no_disturbance_slopes_model 8 1860.6 1882.9 -922.31 1844.6 10.3288 0
## no_interaction_model 11 1861.7 1892.4 -919.86 1839.7 4.9107 3
## full_model_correlated 22 1881.2 1942.5 -918.58 1837.2 2.5592 11
## full_model_uncorrelated 31 1899.3 1985.7 -918.63 1837.3 0.0000 9
## Pr(>Chisq)
## no_system_nr
## no_metaeco_slopes_model 1.0000
## no_disturbance_slopes_model
## no_interaction_model 0.1785
## full_model_correlated 0.9954
## full_model_uncorrelated 1.0000
Odd, in some of the models chi squared is zero and in some the degrees of freedom are zero.
The best model according to AIC should be M + D + (M | day) + (1 |
system number) (no_disturbance_slopes_model). However, it
has zero degrees of freedom.
So, the models I’m going to consider as possible models are the following ones:
Here we want to fit how biomass changes across time to a function. The biomass of our meta-ecosystems looks like this.
To fit these data, we need to produce (and parameterise) a function that resemble how the biomass first increases and then decreases.
Hank produced the following function:
\[biomass = a_4 * (day - a_5) * e^{a_1(day - a_5)}\]
If we parameterise the function and then plot, it looks as follows.
a1 = -0.1
a4 = 1200
a5 = -1
day = seq(0, 30, 0.01)
biomass = a4*(day-a5) * exp(a1*(day-a5))
plot(biomass ~ day)
Now, let’s find the best parameters (a1, a4, a5) that fit our data.
model$m$getPars()
## a1 a4 a5
## -0.1336182 1216.4268967 -1.6147345
And now let’s plot the function to see how it fits our data.
Let’s then include our predictions into the data set.
ds_regional_predicted_shrunk_type = ds_regional %>%
mutate(predicted_from_time = a4*(day-a5)*exp(a1*(day-a5))) %>%
filter(metaecosystem_type == "S_L" | metaecosystem_type == "M_M")
Let’s now work using the fitted data that we found in the previous section.
ds_regional_predicted_shrunk_type_n_day = ds_regional_predicted_shrunk_type %>%
filter(time_point >= 2)
Let’s start from the full model.
model_slopes_correlated = lmer(regional_mean_bioarea ~
predicted_from_time +
metaecosystem_type +
disturbance +
predicted_from_time * metaecosystem_type +
predicted_from_time * disturbance +
metaecosystem_type * disturbance +
predicted_from_time * metaecosystem_type * disturbance +
(predicted_from_time | system_nr),
data = ds_regional_predicted_shrunk_type_n_day,
REML = FALSE)
## Warning: Some predictor variables are on very different scales: consider
## rescaling
(1) Should we keep the correlation in (T | ID)?
model_slopes_uncorrelated = lmer(regional_mean_bioarea ~
predicted_from_time +
metaecosystem_type +
disturbance +
predicted_from_time * metaecosystem_type +
predicted_from_time * disturbance +
metaecosystem_type * disturbance +
predicted_from_time * metaecosystem_type * disturbance +
(predicted_from_time || system_nr),
data = ds_regional_predicted_shrunk_type_n_day,
REML = FALSE,
control = lmerControl(optimizer = "Nelder_Mead"))
## Warning: Some predictor variables are on very different scales: consider
## rescaling
anova(model_slopes_uncorrelated, model_slopes_correlated)
## Data: ds_regional_predicted_shrunk_type_n_day
## Models:
## model_slopes_uncorrelated: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance + predicted_from_time * metaecosystem_type + predicted_from_time * disturbance + metaecosystem_type * disturbance + predicted_from_time * metaecosystem_type * disturbance + ((1 | system_nr) + (0 + predicted_from_time | system_nr))
## model_slopes_correlated: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance + predicted_from_time * metaecosystem_type + predicted_from_time * disturbance + metaecosystem_type * disturbance + predicted_from_time * metaecosystem_type * disturbance + (predicted_from_time | system_nr)
## npar AIC BIC logLik deviance Chisq Df
## model_slopes_uncorrelated 11 1841.8 1872.5 -909.90 1819.8
## model_slopes_correlated 12 1842.7 1876.2 -909.35 1818.7 1.1016 1
## Pr(>Chisq)
## model_slopes_uncorrelated
## model_slopes_correlated 0.2939
No (according to AIC and BIC - however, AIC is only slightly lower).
(2) Should we keep T * M?
model_no_TM = lmer(regional_mean_bioarea ~
predicted_from_time +
metaecosystem_type +
disturbance +
predicted_from_time * disturbance +
metaecosystem_type * disturbance +
predicted_from_time * metaecosystem_type * disturbance +
(predicted_from_time || system_nr),
data = ds_regional_predicted_shrunk_type_n_day,
REML = FALSE,
control = lmerControl (optimizer = "Nelder_Mead")
)
## Warning: Some predictor variables are on very different scales: consider
## rescaling
anova(model_slopes_uncorrelated, model_no_TM)
## Data: ds_regional_predicted_shrunk_type_n_day
## Models:
## model_slopes_uncorrelated: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance + predicted_from_time * metaecosystem_type + predicted_from_time * disturbance + metaecosystem_type * disturbance + predicted_from_time * metaecosystem_type * disturbance + ((1 | system_nr) + (0 + predicted_from_time | system_nr))
## model_no_TM: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance + predicted_from_time * disturbance + metaecosystem_type * disturbance + predicted_from_time * metaecosystem_type * disturbance + ((1 | system_nr) + (0 + predicted_from_time | system_nr))
## npar AIC BIC logLik deviance Chisq Df
## model_slopes_uncorrelated 11 1841.8 1872.5 -909.9 1819.8
## model_no_TM 11 1841.8 1872.5 -909.9 1819.8 0 0
## Pr(>Chisq)
## model_slopes_uncorrelated
## model_no_TM
This is really wired. Df = 0. Let’s then say no then.
(3) Should we keep T * D?
model_no_TD = lmer(regional_mean_bioarea ~
predicted_from_time +
metaecosystem_type +
disturbance +
metaecosystem_type * disturbance +
predicted_from_time * metaecosystem_type * disturbance +
(predicted_from_time || system_nr),
data = ds_regional_predicted_shrunk_type_n_day,
REML = FALSE,
control = lmerControl (optimizer = "Nelder_Mead"))
## Warning: Some predictor variables are on very different scales: consider
## rescaling
anova(model_no_TM, model_no_TD)
## Data: ds_regional_predicted_shrunk_type_n_day
## Models:
## model_no_TM: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance + predicted_from_time * disturbance + metaecosystem_type * disturbance + predicted_from_time * metaecosystem_type * disturbance + ((1 | system_nr) + (0 + predicted_from_time | system_nr))
## model_no_TD: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance + metaecosystem_type * disturbance + predicted_from_time * metaecosystem_type * disturbance + ((1 | system_nr) + (0 + predicted_from_time | system_nr))
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model_no_TM 11 1841.8 1872.5 -909.9 1819.8
## model_no_TD 11 1841.8 1872.5 -909.9 1819.8 0 0
This is really wired. Df = 0. Let’s then say no then.
(4) Should we keep M * D?
model_no_MD = lmer(regional_mean_bioarea ~
predicted_from_time +
metaecosystem_type +
disturbance +
predicted_from_time * metaecosystem_type * disturbance +
(predicted_from_time || system_nr),
data = ds_regional_predicted_shrunk_type_n_day,
REML = FALSE,
control = lmerControl (optimizer = "Nelder_Mead"))
## Warning: Some predictor variables are on very different scales: consider
## rescaling
anova(model_no_TD, model_no_MD)
## Data: ds_regional_predicted_shrunk_type_n_day
## Models:
## model_no_TD: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance + metaecosystem_type * disturbance + predicted_from_time * metaecosystem_type * disturbance + ((1 | system_nr) + (0 + predicted_from_time | system_nr))
## model_no_MD: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance + predicted_from_time * metaecosystem_type * disturbance + ((1 | system_nr) + (0 + predicted_from_time | system_nr))
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model_no_TD 11 1841.8 1872.5 -909.9 1819.8
## model_no_MD 11 1841.8 1872.5 -909.9 1819.8 0 0
This is really wired. Df = 0. Let’s then say no then.
(5) Should we keep T * M * D?
model_no_TMD = lmer(regional_mean_bioarea ~
predicted_from_time +
metaecosystem_type +
disturbance +
(predicted_from_time || system_nr),
data = ds_regional_predicted_shrunk_type_n_day,
REML = FALSE,
control = lmerControl (optimizer = "Nelder_Mead"))
anova(model_no_MD, model_no_TMD)
## Data: ds_regional_predicted_shrunk_type_n_day
## Models:
## model_no_TMD: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance + ((1 | system_nr) + (0 + predicted_from_time | system_nr))
## model_no_MD: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance + predicted_from_time * metaecosystem_type * disturbance + ((1 | system_nr) + (0 + predicted_from_time | system_nr))
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model_no_TMD 7 1844.3 1863.8 -915.15 1830.3
## model_no_MD 11 1841.8 1872.5 -909.90 1819.8 10.49 4 0.03294 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Yes according to AIC (but no according to BIC).
(6) Should we keep (T || ID)?
model_no_random_slope = lmer(regional_mean_bioarea ~
predicted_from_time +
metaecosystem_type +
disturbance +
(1 | system_nr),
data = ds_regional_predicted_shrunk_type_n_day,
REML = FALSE,
control = lmerControl (optimizer = "Nelder_Mead"))
anova(model_no_MD, model_no_random_slope)
## Data: ds_regional_predicted_shrunk_type_n_day
## Models:
## model_no_random_slope: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance + (1 | system_nr)
## model_no_MD: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance + predicted_from_time * metaecosystem_type * disturbance + ((1 | system_nr) + (0 + predicted_from_time | system_nr))
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model_no_random_slope 6 1842.3 1859.0 -915.15 1830.3
## model_no_MD 11 1841.8 1872.5 -909.90 1819.8 10.49 5 0.06249 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Yes according to AIC (but no according to BIC).
(7) Should we introduce the correlation (T | ID) ?
model_no_MD_correlated = lmer(regional_mean_bioarea ~
predicted_from_time +
metaecosystem_type +
disturbance +
predicted_from_time * metaecosystem_type * disturbance +
(predicted_from_time | system_nr),
data = ds_regional_predicted_shrunk_type_n_day,
REML = FALSE,
control = lmerControl (optimizer = "Nelder_Mead"))
## Warning: Some predictor variables are on very different scales: consider
## rescaling
anova(model_no_MD, model_no_MD_correlated)
## Data: ds_regional_predicted_shrunk_type_n_day
## Models:
## model_no_MD: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance + predicted_from_time * metaecosystem_type * disturbance + ((1 | system_nr) + (0 + predicted_from_time | system_nr))
## model_no_MD_correlated: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance + predicted_from_time * metaecosystem_type * disturbance + (predicted_from_time | system_nr)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model_no_MD 11 1841.8 1872.5 -909.90 1819.8
## model_no_MD_correlated 12 1841.5 1875.0 -908.76 1817.5 2.286 1 0.1305
Yes according to AIC (but no according to BIC). This is odd as at the beginning the model told us we shouldn’t keep it.
Should we keep these parameters?
| Parameter | AIC | BIC |
|---|---|---|
| (T | ID) (correlation) | No (but only slightly lower) | No |
| T * M | ? | ? |
| T * D | ? | ? |
| M * D | ? | ? |
| T * M * D | Yes | No |
| (T || ID) | Yes | No |
According to the previous table, then BIC should be the highest for a fixed effect model. Let’s see if this is the case.
fixed_model = lm(regional_mean_bioarea ~
predicted_from_time +
metaecosystem_type +
disturbance,
data = ds_regional_predicted_shrunk_type_n_day)
To make sure everything is okay, let’s take a look at all our models.
anova(model_slopes_correlated, model_slopes_uncorrelated, model_no_random_slope, fixed_model) #Only with models df > 0
## Data: ds_regional_predicted_shrunk_type_n_day
## Models:
## fixed_model: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance
## model_no_random_slope: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance + (1 | system_nr)
## model_slopes_uncorrelated: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance + predicted_from_time * metaecosystem_type + predicted_from_time * disturbance + metaecosystem_type * disturbance + predicted_from_time * metaecosystem_type * disturbance + ((1 | system_nr) + (0 + predicted_from_time | system_nr))
## model_slopes_correlated: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance + predicted_from_time * metaecosystem_type + predicted_from_time * disturbance + metaecosystem_type * disturbance + predicted_from_time * metaecosystem_type * disturbance + (predicted_from_time | system_nr)
## npar AIC BIC logLik deviance Chisq Df
## fixed_model 5 1841.7 1855.6 -915.84 1831.7
## model_no_random_slope 6 1842.3 1859.0 -915.15 1830.3 1.3783 1
## model_slopes_uncorrelated 11 1841.8 1872.5 -909.90 1819.8 10.4897 5
## model_slopes_correlated 12 1842.7 1876.2 -909.35 1818.7 1.1016 1
## Pr(>Chisq)
## fixed_model
## model_no_random_slope 0.24038
## model_slopes_uncorrelated 0.06249 .
## model_slopes_correlated 0.29393
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#anova(model_slopes_correlated, model_slopes_uncorrelated, model_no_TM, model_no_TD, model_no_MD, model_no_TMD, model_no_random_slope, fixed_model) #With all the models
The best model according to AIC is T + M + D + (T || ID)
(model_no_MD), as we found before. However, all other
models seem to be of similar goodness. According to BIC, the best is
(fixed_model)
Therefore, the models I will consider are:
*DOESN’t WORK THOUGH (number of observations (=40) <= number of random effects (=40) )
We want to know if there was a systematic bias in the evaporation of different treatments (disturbance, patch size) and whether evaporation changed across time. My expectation would be that we would see a difference among the exchanges 2,3 and the exchanges 4,5,6. This is because in exchange 2,3 cultures were microwaved in 15 tubes for 3 minutes and in exchange 4,5,6 cultures were microwaved in 4 tubes for only 1 minute.
#Columns: exchange & evaporation
ds_for_evaporation = gather(ds_for_evaporation,
key = exchange,
value = evaporation,
water_add_after_t2:water_add_after_t6)
ds_for_evaporation[ds_for_evaporation == "water_add_after_t2"] = "2"
ds_for_evaporation[ds_for_evaporation == "water_add_after_t3"] = "3"
ds_for_evaporation[ds_for_evaporation == "water_add_after_t4"] = "4"
ds_for_evaporation[ds_for_evaporation == "water_add_after_t5"] = "5"
ds_for_evaporation[ds_for_evaporation == "water_add_after_t6"] = "6"
ds_for_evaporation$evaporation[ds_for_evaporation$exchange == 2] = ds_for_evaporation$evaporation[ds_for_evaporation$exchange == 2] / 2 #This is because exchange contained the topping up of two exchanges
ds_for_evaporation$evaporation[ds_for_evaporation$exchange == 2] = ds_for_evaporation$evaporation[ds_for_evaporation$exchange == 2] + 2 #We need to add 2 ml to the evaporation that happened at the exchange events 1 and 2. This is because we already added 1 ml of water at exchange 1 and 1 ml of water at exchange 2.
#Column: nr_of_tubes_in_rack
ds_for_evaporation$nr_of_tubes_in_rack[ds_for_evaporation$exchange == 1] = 15
ds_for_evaporation$nr_of_tubes_in_rack[ds_for_evaporation$exchange == 2] = 15
ds_for_evaporation$nr_of_tubes_in_rack[ds_for_evaporation$exchange == 3] = 15
ds_for_evaporation$nr_of_tubes_in_rack[ds_for_evaporation$exchange == 4] = 4
ds_for_evaporation$nr_of_tubes_in_rack[ds_for_evaporation$exchange == 5] = 4
ds_for_evaporation$nr_of_tubes_in_rack[ds_for_evaporation$exchange == 6] = 4
It seems like there is no real difference across time, disturbance, or patch type. However, we could also run a mixed effect model to show that they do not.
I’m not going to run it becuase it takes a long time.
mixed.model = lmer(evaporation ~ exchange +
patch_size +
disturbance +
exchange*patch_size +
exchange*disturbance +
patch_size*disturbance +
exchange*patch_size*disturbance +
(exchange||culture_ID) +
(patch_size||culture_ID) +
(disturbance||culture_ID) +
(exchange*patch_size||culture_ID) +
(exchange*disturbance||culture_ID) +
(patch_size*disturbance||culture_ID) +
(exchange*patch_size*disturbance||culture_ID),
data=ds_for_evaporation,
REML = FALSE)
null.model = lm(evaporation ~
1,
data = ds_for_evaporation)
anova(mixed.model, null.model)
To build the mixed effect models we will use the R package lme4. See page 6 of this PDF to know more about the syntaxis of this package.
To do model diagnostics of mixed effect models, I’m going to look at the following two plots (as suggested by Zuur et al. (2009), page 487):
Quantile-quantile plots
Partial residual plots
The effect size of the explaining variables is calculated in the
mixed effect models as marginal and conditional r squared. The marginal
r squared is how much variance is explained by the fixed effects. The
conditional r squared is how much variance is explained by the fixed and
the random effects. The marginal and conditional r squared are
calculated using the package MuMIn. The computation is
based on the methods of Nakagawa, Johnson, and
Schielzeth (2017). For the coding and interpretation of these r
squared check the documentation
for the r.squaredGLMM function
I am starting from the largest model and then simplifying because … (see statistical modelling course at ETH).
I am going to select the best model according to AIC. BIC is better for understanding and AIC for predicting. Halsey (2019) also suggests this approach instead of p values. I’m going to use AIC because I’m interested in knowign how much meta-ecosystem type contributed to the overall regional biomass.
…
culture_info = read.csv(here("data", "PatchSizePilot_culture_info.csv"), header = TRUE)
load(here("data", "morphology", "t0.RData"));t0 = morph_mvt
load(here("data", "morphology", "t1.RData"));t1 = morph_mvt
load(here("data", "morphology", "t2.RData"));t2 = morph_mvt
load(here("data", "morphology", "t3.RData"));t3 = morph_mvt
load(here("data", "morphology", "t4.RData"));t4 = morph_mvt
load(here("data", "morphology", "t5.RData"));t5 = morph_mvt
load(here("data", "morphology", "t6.RData"));t6 = morph_mvt
load(here("data", "morphology", "t7.RData"));t7 = morph_mvt
rm(morph_mvt)
### --- Tidy t0 - t7 data-sets --- ###
#Column: time
t0$time = NA
t1$time = NA
#Column: replicate_video
t0$replicate_video[t0$file == "sample_00001"] = 1
t0$replicate_video[t0$file == "sample_00002"] = 2
t0$replicate_video[t0$file == "sample_00003"] = 3
t0$replicate_video[t0$file == "sample_00004"] = 4
t0$replicate_video[t0$file == "sample_00005"] = 5
t0$replicate_video[t0$file == "sample_00006"] = 6
t0$replicate_video[t0$file == "sample_00007"] = 7
t0$replicate_video[t0$file == "sample_00008"] = 8
t0$replicate_video[t0$file == "sample_00009"] = 9
t0$replicate_video[t0$file == "sample_00010"] = 10
t0$replicate_video[t0$file == "sample_00011"] = 11
t0$replicate_video[t0$file == "sample_00012"] = 12
t1$replicate_video = 1 #In t1 I took only 1 video/culture
t2$replicate_video = 1 #In t2 I took only 1 video/culture
t3$replicate_video = 1 #In t3 I took only 1 video/culture
t4$replicate_video = 1 #In t4 I took only 1 video/culture
t5$replicate_video = 1 #In t5 I took only 1 video/culture
t6 = t6 %>% rename(replicate_video = replicate)
t7 = t7 %>% rename(replicate_video = replicate)
### --- Create ds dataset --- ###
long_t0 = t0 %>% slice(rep(1:n(), max(culture_info$culture_ID)))
ID_vector = NULL
ID_vector_elongating = NULL
for (ID in 1:max(culture_info$culture_ID)){
ID_vector = rep(ID, times = nrow(t0))
ID_vector_elongating = c(ID_vector_elongating, ID_vector)
}
long_t0$culture_ID = ID_vector_elongating
t0 = merge(culture_info,long_t0, by="culture_ID"); rm(long_t0)
t1 = merge(culture_info,t1,by="culture_ID")
t2 = merge(culture_info,t2,by="culture_ID")
t3 = merge(culture_info,t3,by="culture_ID")
t4 = merge(culture_info,t4,by="culture_ID")
t5 = merge(culture_info,t5,by="culture_ID")
t6 = merge(culture_info,t6,by="culture_ID")
t7 = merge(culture_info,t7,by="culture_ID")
ds = rbind(t0, t1, t2, t3, t4, t5, t6, t7); rm(t0, t1, t2, t3, t4, t5, t6, t7)
### --- Tidy ds data-set --- ###
#Column: day
ds$day = ds$time_point;
ds$day[ds$day=="t0"] = "0"
ds$day[ds$day=="t1"] = "4"
ds$day[ds$day=="t2"] = "8"
ds$day[ds$day=="t3"] = "12"
ds$day[ds$day=="t4"] = "16"
ds$day[ds$day=="t5"] = "20"
ds$day[ds$day=="t6"] = "24"
ds$day[ds$day=="t7"] = "28"
ds$day = as.numeric(ds$day)
#Column: time point
ds$time_point[ds$time_point=="t0"] = 0
ds$time_point[ds$time_point=="t1"] = 1
ds$time_point[ds$time_point=="t2"] = 2
ds$time_point[ds$time_point=="t3"] = 3
ds$time_point[ds$time_point=="t4"] = 4
ds$time_point[ds$time_point=="t5"] = 5
ds$time_point[ds$time_point=="t6"] = 6
ds$time_point[ds$time_point=="t7"] = 7
ds$time_point = as.character(ds$time_point)
#Column: eco_metaeco_type
ds$eco_metaeco_type = factor(ds$eco_metaeco_type,
levels=c('S', 'S (S_S)', 'S (S_L)', 'M', 'M (M_M)', 'L', 'L (L_L)', 'L (S_L)'))
#Select useful columns
ds = ds %>%
select(culture_ID,
patch_size,
disturbance,
metaecosystem_type,
mean_area,
replicate_video,
day,
metaecosystem,
system_nr,
eco_metaeco_type)
#Reorder columns
ds = ds[, c("culture_ID",
"system_nr",
"disturbance",
"day",
"patch_size",
"metaecosystem",
"metaecosystem_type",
"eco_metaeco_type",
"replicate_video",
"mean_area")]
As in Jacquet, Gounand, and Altermatt (2020) I will create 12 size classes.
#### --- PARAMETERS & INITIALISATION --- ###
nr_of_size_classes = 12
largest_size = max(ds$mean_area)
size_class_width = largest_size/nr_of_size_classes
size_class = NULL
### --- CREATE DATASET --- ###
size_class_boundaries = seq(0, largest_size, by = size_class_width)
for (class in 1:nr_of_size_classes){
bin_lower_limit = size_class_boundaries[class]
bin_upper_limit = size_class_boundaries[class+1]
size_input = (size_class_boundaries[class] + size_class_boundaries[class + 1])/2
size_class[[class]] = ds%>%
filter(bin_lower_limit <= mean_area) %>%
filter(mean_area <= bin_upper_limit) %>%
group_by(culture_ID,
system_nr,
disturbance,
day,
patch_size,
metaecosystem,
metaecosystem_type,
eco_metaeco_type,
replicate_video) %>% #Group by video
summarise(mean_abundance_across_videos = n()) %>%
group_by(culture_ID,
system_nr,
disturbance,
day,
patch_size,
metaecosystem,
metaecosystem_type,
eco_metaeco_type) %>% #Group by ID
summarise(abundance = mean(mean_abundance_across_videos)) %>%
mutate(log_abundance = log(abundance)) %>%
mutate(size_class = class) %>%
mutate(size = size_input) %>%
mutate(log_size = log(size))
}
ds_classes = rbind(size_class[[1]], size_class[[2]], size_class[[3]], size_class[[4]],
size_class[[5]], size_class[[6]], size_class[[7]], size_class[[8]],
size_class[[9]], size_class[[10]], size_class[[11]], size_class[[12]],)
The code for plotting the size class distribution across small patches doesn’t work for some reason. It displays all plots as if they were the same.
## [1] "Disturbance = low"
## [1] "Disturbance = high"
Regional bioarea density (mean bioarea density between two patches) in meta-ecosystems of the same total area, but whose two patches have either the same size (red, medium-medium meta-ecosystem) or that have a smaller and larger patch (blue, small-large meta-ecosystem). Points represent the mean, error bars represent the standard deviation. For clarity, only the low disturbance treatment is shown here. See the Appendix for equivalent the figure of the high disturbance treatment.
Local biomass production in patches that are either connnected to a patch of the same size (green, small patches of small-small meta-ecosystemes) or to a patch of larger size (orange, small patches of small-large meta-ecosystems). Points represent the mean, error bars represent the standard deviation. For clarity, only the low disturbance treatment is shown here. See the Appendix for equivalent the figure of the high disturbance treatment. Problem: the small patches in green actually have different sample size than the one in orange. Does it still make sense to include the standard deviation?
## Time difference of 29.68097 secs